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Abstract 

We introduce a new class of singular partial differential equations, referred to as the 
second-order hyperbolic Fuchsian systems, and we investigate the associated initial value 
problem when data are imposed on the singularity. First, we establish a general existence 
theory of solutions with asymptotic behavior prescribed on the singularity, which relies on 
a new approximation scheme, suitable also for numerical purposes. Second, this theory 
I is applied to the (vacuum) Einstein equations for Gowdy spacetimes, and allows us to 

. recover, by more direct arguments, well-posedness results established earlier by Rendall 

and collaborators. Another main contribution in this paper is the proposed approximation 
scheme, which we refer to as the Fuchsian numerical algorithm and is shown to provide 
i | highly accurate numerical approximations to the singular initial value problem. For the 

class of Gowdy spacetimes, the numerical experiments presented here show the interest 
and efficiency of the proposed method and demonstrate the existence of a class of Gowdy 
5^ spacetimes containing a smooth, incomplete, and non-compact Cauchy horizon. 
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1 Introduction 

In general relativity and, more generally, in the theory of partial differential equations, singular 
solutions play a central role in driving the development of the theory and interpretation of the 
results - from the discovery of the Schwarzschild solution to the Penrose and Hawking's celebrated 
singularity theorems (cf. for instance [24]), as well as in the modern efforts to understand the 
cosmic censorship and BKL (Belinsky-Khalatnikov-Lifshitz) conjectures [H HU] . 

In fundamental and pioneer work, Choquet-Bruhat |10j established that the initial value 
problem associated with Einstein's field equations is well-posed (in suitable Sobolev spaces) 
and, later together with Geroch [11], that for each choice of initial data consistent with the 
constraints there exists a unique maximal globally hyperbolic development. This theory allows 
to define a unique correspondence, between the space of solutions of the Einstein equations and 
the space of initial data set on a given (spacelike) hypersurface. In principle, solutions could 
be extended until a singularity forms and this theory provides the basis to tackle outstanding 
questions about the long-time behavior of solutions, in both future and past directions. Consider 
for instance the situation that singularities develop in the past of the given hypersurface, and let 
us refer to the construction based on [TO] |TT| as the "backward approach" . In practice, however, 
revealing information about singular solutions requires additional techniques of analysis going 
much beyond [101 lllj . On one hand, there is no a-priori information whether a given choice of 
initial data evolves into a "singularity" at all. On the other hand, if a singularity arises, it is not 
a-priori clear at "which location" it will occur, nor what kind of "singular solution" it will be. 
This approach requires a global-in-time control of solutions which is hopeless in many cases due 
to the complexity of the nonlinear field equations, the freedom of choice of gauge etc. In fact, 
the "backward approach" has so far led to successful conclusions only under strong symmetry 
assumptions; cf. Ringstrom [30] (and the references therein). 

The so-called Fuchsian method, which we refer to here as a "forward" approach, provides 
an alternative to the above (backward) approach. It has the advantage of being simpler to deal 
with in practice (both analytically and numerically) in the class of problems it does apply, but 
has the disadvantage that it does not allow - for essential reasons - to handle the full solution 
space and the conclusions are stated in terms of "generic data" . Nevertheless it allows to study 
particular classes of singular solutions. The main idea is to give data "on the singularity" , that 
is, to prescribe the leading-order behavior of solutions in a neighborhood of the singular time 
and, then, to evolve forward in time away from the singularity - in contrast to prescribing data 
on a Cauchy surface and evolving towards the singularity. (We will make this idea precise below.) 
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Our main objective in the present paper is precisely to further study this singular initial value 
problem (corresponding to the forward direction). 

The Fuchsian method was introduced to general relativity by Kichenassamy and Rendall [28] 
who covered a class of singular equations, the so-called "Fuchsian partial differential equations" , 
while establishing that Einstein's field equations under Gowdy symmetry are included in this 
class. Under suitable analyticity conditions on the solutions they proved the well-posedness 
of the singular initial value problem (without imposing any hyperbolicity condition). Later, 
Rendall |36j took the hyperbolicity property into account and generalized the theory to the 
smooth solutions to the Gowdy equations. For further results about Fuchsian equations we refer 
to [2UI H"2l [TBI IT4"] ; especially, in [13], a generalization of the standard Fuchsian theory was recently 
introduced. 

Our motivation in this paper now is threefold. First, it is of general interest to find a 
reliable numerical scheme for the singular initial value problem of hyperbolic Fuchsian equations, 
which would allow us to construct solutions to Einstein's field equations with prescribed singular 
behavior. This was indeed the motivation of Amorim, Bernardi, and LeFloch in [T]. However, the 
approximation scheme suggested by the theoretical works above is neither suitable nor natural 
for the numerical treatment, and it is not easy to obtain error estimates which are necessary in 
practice to judge the quality of the numerical solutions. In this paper, we address this drawback 
by introducing a new approximation scheme. 

Secondly, we want to perform these studies — both theoretically and numerically — as eco- 
nomically as possible. Importantly, the Einstein equations for Gowdy spacetimes are naturally 
expressed in a second-order form. This motivates us to develop here a theory of second-order 
Fuchsian equations which is advantageous as it saves us from the hassle of turning the system 
into a first-order form first, as required by the classical Fuchsian theory. In fact, there is no unique 
way of turning a second-order system into a first-order form, and this issue can be particularly 
problematic when the solutions are expected to be singular. It has also been recognized in the 
numerical literature that the direct discretization in second-order form leads to more accurate 
results [3"T] . 

Our third motivation for this paper is that the existing works make no statement about how 
to guess the leading-order part. This is a delicate issue and we prove here that, for a large class 
of systems and in way compatible with the (already mentioned) BKL conjecture, it is possible 
to make a "canonical guess" , as we explain below. 

These three main issues are addressed in the present paper, which is organized as follows. 

Sections [2] and [3] are the main theoretical part. In Section |2.1[ we introduce the class of 
equations of interest in this work. We focus on a class of hyperbolic Fuchsian equations, the 
second-order hyperbolic Fuchsian equations, which are systems of semi-linear wave equa- 
tions including certain singular terms. This section includes a rigorous definition of the singular 
initial value problem for this class of equations. Section 2J5 is devoted to determine the leading- 
order behavior of solutions and hence to make a "canonical guess". We study the case when 
the principal part of the equation dominates over the source-term at t = in a certain sense 
and derive a canonical two-term expansion. The reasoning is based on suitable heuristics 
compatible with the BKL conjecture. (Later in this paper we derive precise conditions under 
which this heuristics is justified.) 

The singular initial value problem is discussed rigorously in Section [3] We introduce our new 
approximation scheme which, both, yields a simple and direct proof of existence of solutions to 
the singular initial value problem, and a natural numerical algorithm (introduced in Section 3.4) 
for which practical error estimates can be obtained. The main idea is to approximate the solution 
of the singular initial value problem by a sequence of solutions to the standard (regular) initial 
value problem. 
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Finally, in Section [4j we turn our attention to the class of Gowdy spacetimes satisfying 
Einstein's field equations which is an application of particular interest, and we demonstrate the 
practical use of our theory. First, we recapitulate the standard heuristic arguments for Gowdy 
spacetimes and demonstrate that these are consistent with the heuristics introduced earlier (with 
some important subtleties discussed below). In Section 4.3 the now classical results by Rendall 
[3"rj] are recovered while our new approach shed some new light on the Gowdy equations. In 
Section [5] we present numerical solutions to the singular initial value problem associated with the 
Gowdy equations. After some test cases, we numerically construct solutions with incomplete, 
non-compact Cauchy horizons. 



2 Second-order hyperbolic Fuchsian systems 
2.1 A class of singular equations 

Definition 2.1. A second-order hyperbolic Fuchsian system is a set of partial differential equa- 
tions of the form 

D 2 v + 2ADv + Bv-t 2 K 2 d 2 x v = f[v], (2.1) 

in which the function v : (0, 5] x U — > R" is the main unknown (defined for some S > and some 
interval U ), while the coefficients A = A(x), B = B(x), K — K(t, x) are diagonal n x n matrix- 
valued maps and are smooth in x £ U and t in the half-open interval (0,5], and f = f[v](t,x) is 
an n-vector-valued map of the following form 

f[v](t, x) := f(t, x, v(t, x), Dv(t, x), tK(t, x)d x v(t, x)j . 

We assume that the time variable t satisfies t > and use the operator D := tdt to write the 
equations. The equation is henceforth assumed to be singular at t — 0. The assumption that U 
is a one-dimensional domain makes the presentation simpler, but most results below remain valid 
for arbitrary spatial dimensions. For definiteness and without much loss of generality, we assume 
throughout this paper that all functions under consideration are periodic in the spatial variable x 
and that U is the periodicity domain. All data and solutions are extended by periodicity outside 
the interval U. Moreover, we assume that the coefficients A and B do not depend on t, see below. 
We denote the eigenvalues of A and B by . . . , a^ and . . . , b^ n \ respectively. When it 
is not necessary to specify the superscripts, we just write a, b to denote any eigenvalues of A, B. 
With this convention, we introduce: 

Ai := a + \Jd 2 - b, A 2 := a - y/a 2 -b. (2.2) 
It will turn out that these coefficients, which might be complex in general, are important to 



describe the expected behavior at t = of general solutions to (2.1 1. Further restrictions on 
the coefficients and on the right-hand side will be imposed and discussed in the course of our 
investigation. 

After a suitable reduction to a first-order form, our choice of singular equations falls into 
the class of hyperbolic Fuchsian equations 28, 3|3 [T3]. We make this particular choice of 
equations here for the following reasons. First we restrict to hyperbolic equations because this 
is the case of interest and allows us to control solutions in much greater detail than without this 
assumption. Second, the equations are kept in second-order form here because it is economic 
and efficient to do this for the applications we have in mind both analytically and numerically. 
Third, the particular class of equations allows us to simplify the presentation of the general 
results obtained in this paper. However, all results presented here can be generalized to a 
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general class of symmetric hyperbolic Fuchsian equations - compare to the discussion in Rendall 
[36] - for arbitrary spatial dimensions. Even beyond this it is possible to generalize the theory to 
equations with time-dependent coefficients A and B. Also the restriction to spatial periodicity 
is not essential because these equations obey the domain of dependence property just as usual 
non-singular hyperbolic equations under suitable assumptions on K, see below. 

The eigenvalues of the matrix K are denoted by k^> and, in the scalar case (or when there is 
no need to specify the index), we simply write k. These quantities are interpreted as characteristic 
speeds. Throughout this section, we assume that they have the form 

k^(t,x)=t^^^(t,x), 

with /3 (l) : U (-1, oo), v (i) : [0, S] x U -> (0, oo) smooth functions. 

In particular, we assume that each derivative of i>w has a unique finite limit at t = for each 
x E U. Note that we allow for the characteristic speeds to diverge at t = 0. At a first glance, 
this appears to conflict with the standard finite domain of dependence property of hyperbolic 
equations. A closer look at the requirement (3{x) > — 1, however, indicates that the characteristic 
curves are integrable at t = and, hence that the finite domain of dependence property is 
preserved under our assumptions. 

The operator associated with the principal part of the system is 

L := D 2 +2AD + B -t 2 K 2 d 2 x -. L-t 2 K 2 d 2 x . (2.4) 



This is a linear wave operator for t > and, indeed, (2.1 ) is hyperbolic for all t > 0. Later on we 
will construct solutions where the first three terms L of the principal part are of the same order 
at t = and "dominant" , while the source-term of the equation as well as the second spatial 
derivative term are assumed to be of higher-order in t at t = and hence "negligible" . Note that 
at this level generality, there is some freedom in bringing terms from the principal part to the 
right-hand side of the equation, and absorbing them into the source- function / (or vice- versa). 
This freedom has several (interesting) consequences: roughly speaking, some normalization will 
be necessary later, yet at this stage, we do not fix the behavior of / at t = 0. 

2.2 Singular initial value problem 

Consider any second-order hyperbolic Fuchsian system with coefficients a, b, Ai,A2, satisfying 



(2.2). To simplify the presentation, we restrict attention to scalar equations (n = 1) and shortly 
comment on the general case in the course of the discussion. 

Fix some integers l,m > and constants a, S > 0. For w £ C l ((0, 5],H m (U)), we define the 
norm 

i , 1/2 



W 



su p \yy>i / t mx2(x) ~ a) \d^D p W (t,x)\ 2 dx 

°<*<* Vp=0g=0^ 



and denote by Xs t a,i,m the space of all such functions with finite norm || w\\s tC t,l,m < oo. Through- 
out, H m (U) denotes the standard Sobolev space and we recall that all functions are periodic 
in the variable x with U being a periodicity domain. To cover a system of n > 1 second-order 
Fuchsian equations, the norm above is defined by summing over all vector components with 
different exponents used for different components. Recall that each equation in the system will 
have a different root function A 2 . We allow that a — (a^\ . . . ,a^ n ') is a vector of different 
positive constants for each equation. The constant 5, however, is assumed to be common for all 
equations in the system. With this modification, all results in the present section remain valid 
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for systems of equations. We comment later that in fact, a is not required to be a constant, but 
in most of the following results it will be treated like a constant for simplicity. The reason to 
include the quantity A2 into the definition of the norms is motivated by the canonical choice of 
the leading-order term introduced later. Throughout it is assumed that 5RA2 is continuous and 
it is then easy to check that {Xs, a ,l,mi II ■ \\s,a,i,m) is a Banach space. 

For the discussion of hyperbolic equations, it makes sense to also introduce the following 
Banach spaces. For each non-negative integer / and real numbers 5, a > 0, we define Xs, a ,l '■— 
f1p=o Xs,a,p,i-p : and introduce the norm 



1/2 

l/llr5,a,p,i-p I 1 / ^-S.a, 



\p=0 



As we will see in the course of the following, however, it is not possible to control solutions of 
our equations in the spaces Xg <a i directly. It turns out that we must use spaces (X$ a j, \\ ■ Wg 1 ; ) 
instead. These are defined as earlier, but in the norm ||/|| ^ Q ; of some function /, the highest 
spatial derivative term d l x f is weighted with the additional factor t@ +1 . Here (3 is the exponent 



of the characteristic speed given by (2.3 1. It is easy to sec under the earlier conditions that 
also (Xs, a ,i, || • Ws'a i) are Banach spaces. We also note that Xs, a ,i C Xs, a ,i- Let us also define 
Xs : a : oo ■— f)i=o-^s,a,h an d note that X StCCi00 = C\ilo -^s,a,i- 
For w £ Xg a i we set 



w 



/OO />00 . -I 

/ t x ^w(s,y)kJlog -)k v (x - y)-dsdy. 
-00 Jo \ t/ s 



Here, : M — > M + is a smooth kernel supported in [—77,77], satisfying J R k v (x)dx = 1 for all 
positive 77. Then w v is an element of Xs y a-e,i H C°°((0, 6} x U) for every e > 0. Furthermore, the 

sequence of such mollified functions in the limit 77 — > converges to w in the norm || • ||^ ; . 

Hence any element in Xg, a ,l can be approximated by smooth functions. 

We mentioned earlier that we are interested in solving the "forward problem" , referred to as 
the singular initial value problem (SIVP) - in this paper. More precisely, we need to guess 



a leading-order term u of solutions v to (2.1 1 so that the remainder 



w(t, x) := v(t, x) — u(t, x), 

can be interpreted as "higher order" in t at t = 0. By this we mean that w is an element in 
Xs,a,i for some (sufficiently large) a > on a small time interval (0,<5]. If for a given u such a 
solution v exists then we say that v obeys the leading-order behavior given by u. Often u 
will be parametrized by certain free functions which we call asymptotic data, see below. For 
later convenience, we introduce the operator F as 

F[w]{t,x):=f[u + w](t,x). (2.5) 



2.3 Canonical leading-order term 

The main first step for solving the singular initial value problem is to guess a leading-order 
term u. In some applications this can be very tricky, but in many situations, which we will be 
most interested in in this paper, one can make a canonical guess. These situations are described 
heuristically as follows. 
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Canonical two-term expansion 



Consider the principal part operator L in (2.4) and note that it incorporates certain lower 



derivative terms. The reason for writing L like this is that we expect in many cases that these 
terms are significant and of leading-order at the singularity t — 0. In contrast, the source- 
term and spatial derivatives can often be anticipated as negligible in some sense under suitable 
assumptions below. This is motivated by the BKL conjecture in general relativity. In order to 
make this more concrete, let us assume that the leading-order term is an exact solution of the 
system of ordinary differential equations (parametrized by x), which is obtained when all terms 
in the equation, except for those given by L, are set to zero. We refer to this leading-order term 
u as the "canonical leading-order term" or the canonical two-term expansion: 

log* + u**(x) *"<*(*), (a(x)f = b(x), 
H,X) \u*(x)t- x ^+u t *(x)t- x ^, (a(x)f^b(x), [ ' 

for some freely prescribed asymptotic data and it** £ H m (U), where m' is some non- 
negative integer. We refer to this as the "Fuchsian heuristics" because the leading-order behavior 
will be determined by Fuchsian ordinary differential equations. 

We clearly see the dependence of the expected leading-order behavior at t = on the coef- 
ficients of the principal part of the equation. If the roots Ai and A2 are real and distinct, i.e. if 
a 2 > b, we expect a power-law behavior. In the degenerate case Ai = A2, i.e. if a 2 = 6, we expect 
a logarithmic behavior. Finally, when Ai and A2 are complex for a 2 < 6, the solution is expected 
to have an oscillatory behavior at t = of the form 

u(t,x) = t~ a( - x \u* cos(Xi(x) logt) + m** sin(A/(x) logt)) + . . . 

for some real coefficient functions u#(x) and u**(x); note that in this case, Xi = A 2 = a + iXj 
with A/ := \Jb 2 — a. 

If the coefficients of the equations are such that there is a continuous transition between 



the two cases in (2.6), then the asymptotic data functions it* and u** must be renormalized as 
follows. Define T(x) := \J a(x) 2 — b(x) which might be real or imaginary dependent on the values 
of the coefficients. If there are points xo £ U so that T(xo) = and other points x\ £ U with 
r(x 1 ) 7^ 0, then let us set 

, , _ M x ) -u^{x)/T(x) _ u*(x)+u**(x)/T(x) 

u^yx) — ^ ' u**\x) — ^ i \ Z -'J 

and choose u**(a:) as asymptotic data functions. This guarantees that u(t, x) given by 



(2.6) is smooth for all t > 0. 



Higher-order canonical expansions 

For some applications we require expansions of the solutions at t = with more than two terms 
in order to describe the leading-order behavior. A particular important example is the Gowdy 
case in Section [4j Following Rendall [35], those can be constructed as follows, without going into 
the details. Consider the Fuchsian ODE case of (12 .11), written here for the scalar case only, 



D 2 v(t,x) + 2a(x)Dv(t, x) + b(x)v(t,x) = f[v](t,x), 

where x is interpreted as a parameter. Let first f[v](t,x) — fo(t,x) be a given function. Under 
suitable decay assumption^ on fa at t — 0, there exists a unique solution v of this equation 

1 Details can be found in [71 18). 
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obeying the canonical two-term expansion u given by (2.6) for given asymptotic data functions 



u* and it**. Let H be the operator mapping /o to the remainder w = v — u of the solution v. 



Now consider an arbitrary source-term f[v] and let the operator F be as defined in (2.5) for the 
given function u. Let w\ = and 

w j+ i := H o F[wj], j G N. 

Finally, set v 3 ■ — u + Wj for all j G N. Clearly, v\ = u. One finds that the order of u,+i — Vj in t 
at t = increases with j. Hence vj can be interpreted as an expansion of the solution at t = 
whose order in t increases with j. Moreover, it turns out that the order of the residual in t at 
t = 0, obtained when Vj is plugged into the equation, increases with j. Thus Vj can be considered 
as an asymptotic solution of the Fuchsian ODE. In many situations it is thus meaningful to use 
Vj as the canonical leading-order term u for any given j G N. 

Limitations of this heuristics 

At this stage it is of course highly unclear under which conditions there exist solutions of the 



equations which obey the leading-order term u given by (2.6) or any of its higher-order version 
Vj constructed before. The resolution of this problem will be central to this paper. In many 
applications in general relativity, the canonical two-term expansion is the correct guess for the 
leading-order term, if the asymptotic data are consistent with constraint equations implied by 
Einstein's field equations. However, we know of several cases when the operator L does not give 
rise to the dominant term at t — 0. For example for the Gowdy case, nonlinear terms from the 
source-term need to be taken into account. As we discuss there, however, the problem can be 
reduced to the canonical case by adding a certain term to the equation. For Gowdy solutions with 
spikes [3 [371 123 131 HI GH3 HD GUI HD] the situation is significantly more complicated because then, 
other nonlinear terms and spatial derivative terms can become significant. Another important 
example is given by the mixmaster dynamics HI [251 Ed SO] ■ There, one has to control 
a complicated interplay between nonlinear terms in the source-term in order to fix the leading- 
order term. In general when the equations are generalized to time-dependent coefficients in the 
principal part or to the quasi-linear case, the notion of canonical two-term expansions will apply 
only under suitable conditions. 

3 Well-posedness theory and the Fuchsian numerical algo- 
rithm 

3.1 An approximation scheme 



We begin with some notation. For w e Xs t a,i, the operator L in (2.4) is defined in the sense of 
distributions, only, via: 



(C[w],(, 



t n\ a (x)-af -Dw(t,x)D</>(t,x) + (2A{x) -5RA 2 (x) +a - 1) Dw(t,x)4>(t,x) 
+ B{x) w(t, x)4>(t, x) + tK(t, x)d x w(t, x)tK(t, x)d x <f)(t, x) 
+ (2td x K(t,x) + d x ^\ 2 {x)K{t 1 x)t\ogt)tK(t 1 x)d x w{t,x) (j)(t,x) \ dxdt, 



where <fi is any test function, i.e. a real-valued C°°-function on (0, S] x R together with some 
T G (0, S) and a compact (i.e. closed and bounded) set K G K so that <j>(t, x) = for all t > T 
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and x K, and each derivative of <fi has a finite (not necessarily vanishing) limit at t = for 
every x € U. For our later discussion, we note that for any given test function (j>, the linear 
functional (£[■], (^) : Xs. a ,i — > K is continuous with respect to the norm || • ||^ Q 1 . This is the 
main reason to include the factor t^^^-a m definition of C. 



If the operator F defined by (2.5) for a given leading-order term u gives rise to a map 
Xs, a ,i Xs^ a ,o, where w H ► .F[iu], it is meaningful to define its weak form by (for all test 
functions <f>) 

(T[w],<j>) '.= f [ t 3iX2 ^- a F[w](t,x)(t)(t,x)dxdt. 
Jo Jr 

Definition 3.1 (Weak solutions of second-order hyperbolic Fuchsian systems). Let u be a given 
function and 8, a > be constants. Then, one says that w € Xs. a ,i is a weak solution to the 



second-order hyperbolic Fuchsian equation (2.1), provided 



V[w] := C[w] + C[u] - F[w] = 0. 

Let us now start our discussion with the linear case of second-order hyperbolic Fuchsian 
equations and introduce our new approximation scheme. The following conditions are assumed: 

1. Vanishing leading-order part: u = 0. 

2. Linear source-term: 

F[w](t,x) = f Q (t,x) + f 1 (t,x)w + f 2 (t,x)Dw + f 3 (t,x)tkd x w, (3.1) 

with given functions /o, /i, /3, so that / 1; / 2 , /3 are smooth spatially periodic on 
(0, 6] x U, and near t = 



sap f a (t,x) = 0(lP), a = 1,2, 3, (3.2) 



for some constant /i > 0. 



We have not made any assumptions for the function /o yet, since in the following discussion this 
function will play a different role than /i, f%, fa. Moreover, no loss of generality is implied by 
the condition u = 0, since the general case can be recovered by absorbing L[u] into the function 
fa- 

Under these assumptions, we pose the question whether there exists a unique weak solution 
w in Xs, a ,i for some S, a > of the given second-order hyperbolic equation. The main idea here 
is our new approximation scheme. We approximate a solution of the singular value problem by 
a sequence of solutions of the regular initial value problem. 

Definition 3.2 (Regular initial value problem (RIVP)). Fix to € (0, 8] and some smooth periodic 



functions g, h : U — > R, and suppose that the right-hand side is of the form (3.1 ) with given smooth 



spatially periodic functions fo, fx, fi, fa on [to, 8] x U. Then, w : [to, 8] X U — > K. is called a 



solution of the regular initial value problem associated with the "regular data" g,h, if (2.1) 
holds everywhere on (to, 8] x U and, moreover, the remainder w :— v — u, for some function u, 
satisfies 

w(t , x) = g(x), d t w(t , x) = h(x). 

For the regular initial value problem, we indeed assume that fo is smooth, just as f\, f-x and 
fa- By the general theory of linear hyperbolic equations, the regular initial value problem is 
well-posed, in the sense that there exists a unique smooth solution w defined on [to, 8] for any 
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choice of smooth regular data. Let u be a choice of leading-order term. Let (t n ) be a sequence 
of positive times converging to zero so that each t n is smaller than some 5 > 0. For each n G N, 
we construct an approximate solutions w n of the singular initial value problem as follows. Let 
w n e on the time interval (0, t n ]. On the time interval [t n ,5] we set w n to be the solution 
of the regular initial value problem with initial time t n and zero regular data. Hence, w n is in 

c 1 ((q,S\xU)nx s , a ,i. 

The central result of this section is that this sequence of approximate solutions converges to 
the solution of the singular initial value problem under suitable conditions. 

Proposition 3.3 (Existence of solutions of the linear singular initial value problem in Xs. a ,i)- 



Under the assumptions (3.1 1 and (3.2 1 7 the sequence (w n ) of approximate solutions with initial 
times (t n ) — > converges to the unique solution w G Xs. a ,\ of the singular initial value problem 
for given 6, a > provided: 

1. The matrix 

{3t(\ 1 -\ 2 ) + a ((SA 1 ) 2 /t 7 -77)/2 

N := {{^Xxf/rj - T))/2 a td x k - djt^ - A 2 )(ifclogt) 

\ td x k - dJR{\i - \ 2 ){tk\ogt) 9ft(Ai - A 2 ) + a - 1 - Dk/k 

(3.3) 

is positive semidefinite at each (t,x) G (0,<5) x U for a constant r\ > 0. 

2. The source-term function /o is in Xs, a +e,a f 0T some e > 0. 

Then, the solution operator H : Xs ta +e,o ~^ Xg iCti i, /q i— > w, is continuous and there exists a 
finite constant C' e > so that 

||H[/ ]||^ )1 <^C' £ ||/o|| 5 , Q+£ ,o ! 

for all such fa. The constant C e can depend on 5, but is bounded for all small S. The approximate 
solutions w n satisfy the following error estimate for all n,m G N: 

\\w m - tu„||£ Q|1 < C\G(t n ) - G(t m )\, 

where C > is a constant and 

G(t) := f s^Ws^-^fois^^L^u)^. 
Jo 

We call N the energy dissipation matrix. We have assumed that a is a positive constant. 
If, however, a is a positive spatially periodic function in C l (U), the definition of the spaces 
Xs.a.k and Xs^ a ,k remains the same, and only the (2, 3)- and (3, 2)-components of the energy 
dissipation matrix TV change to td x k — 9 x (3?(Ai — A 2 ) + a)(tk log t). In the following, we continue 
to assume that a is a constant in order to keep the presentation as simple as possible, but we 
stress that all following results hold (with this slight change of N) if a is a function, and hence 
no new difficulties arise. 

The proof is based on controlling the energy of the approximate solutions. Let us restrict 
our presentation here to the scalar case n — 1 for this whole section; the general case can be 
obtained with the same ideas. Choose S, a > and let w G C 1 ((0, S] X U) be a spatially periodic 
function. Then, we define its energy at the time t G (0, 8] by 



E[w]{t) :=e- Kt ~' / t 2 ^^' a ^e[w}(t,x)dx, 

Ju (3.4) 
e[w](t,x) ((r/w(t,x)) 2 + (Dw(t,x)) 2 + {tk{t,x)d x w(t,x)) 2 ) , 
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for some constants k > 0, 7 > and 77 > 0. For convenience, we also introduce the following 
notation. For any scalar- valued function w, we define the vector- valued function 



w(t, x) := t UX2ix) - a {r]w(t, x),Dw(t, x),tk(t, x)d x w(t, x)), 
involving the same constants as in the energy. Then, we can write 



(3.5) 



E[w](t) = l -e 



\w{t,-)\ 



the norm here being the Euclidean L 2 -norm for vector-valued functions in x. It is important 
to realize that, provided 77 > 0, the expression sup 0<t <^ Oil f° r functions of the form 
(3.5) yields a norm which is equivalent to || • ||^ a 1; thanks to (2.3). Therefore, the energy (3.4) 



is of relevance for the space Xg a ,i- 

Of central importance for the proof of Proposition |3.3| are energy estimates for the regular 
initial value problem. 

Lemma 3.4 (Fundamental energy estimate for the regular initial value problem). Suppose that 



the source term is of the form (3.1) with the condition (3.2) and that the energy dissipation 



matrix (3.3) is positive semidefinite on (0,5] x U for given constants a,r/ > and a sufficiently 
small 5 > 0. Then there exist constants C,k,j > 0, independent of the choice of to € (0,(5], so 
that for all solutions w of the regular initial value problem with smooth regular data at t = t , 
we have 

\\w(t, -)\\ LHU) < Ce§*(* T -'S) (\\w( to , -)\\ LHU) + f s-V^Ms, Wvwds) , 
for all t £ [to, 5]. 

The proof of this lemma is standard. However, one has to confirm that the energy stays 



uniformly finite at t = 0, and this is guaranteed by the positivity condition for N in (3.3). 
Moreover, this results demonstrates the importance of the assumption j3{x) > —1 in (2.3). 
Namely, if /3(x) < — 1 at a point x E U, then for any choice of a and 77, the matrix N would 
not be positive semidefinite for small t at x. While the energy estimate would still be true for a 
given to, we would nevertheless lose uniformity of the constants in the estimates with respect to 
to- We stress that this uniformity is crucial in the proof of Proposition 3.3 All proofs and more 
details can be found in [?]• 



3.2 Nonlinear theory 

In turns out that the space Xs. a ,i is too large for our nonlinear theory. Namely, we will need 
to require a Lipschitz property of the source-term for which this space rules out natural nonlin- 
earities, for instance quadratic ones. Moreover, the statement that the solution of the Fuchsian 
equation w is an element of Xg Q l yields only weak information about the behavior of the first 
spatial derivative at t = 0, which is indeed not sufficient to interpret w and all its first derivatives 
as the remainder of the solution. We resolve these problems by going to Xg^ a k for some k > 1. 
The first issue above disappears for k = 2 in one spatial dimension. In some applications, the 
spaces Xg :a ^k still impose a too strong restriction due to the weak control of the highest spatial 
derivative. In such a case we are required to formulate the theory in the space Xg^ a>ao . 

The first step is to reconsider the linear case and derive a result analogous to Proposition |3.3| 
in the spaces Xg^ a ,k for arbitrary fc by making suitable stronger assumptions on fo. It turns out 
that higher spatial derivatives of the solutions may involve additional logarithmic terms in t. As 



11 



a consequence we find that the energy dissipation matrix must be assumed to be positive definite 
instead of positive semidefinite. With this at hand, the central result of this section, which we 
write for k = 2 for definiteness, is the following. 

Proposition 3.5 (Existence of solutions of the nonlinear singular initial value problem in Xs, a ,2)- 



Suppose that we can choose a > so that the energy dissipation matrix (3.3 1 is positive definite 
at each (t, x) £ (0, 5) x U for a constant n > 0. Suppose that u = and that the operator F 
has the following Lipschitz continuity property: For a constant e > and all sufficiently small 
5, the operator F maps Xs -at 2 into Xs ia +e,i and, moreover, for each r > there exists C > 
( independent of 5) so that 

\\F[w] - F[w]\\ Sta+e>1 <C\\w- w\\s <at2 (3.6) 



for all w,w £ B r (0) C Xs, a ,2- Then, there exists a unique solution w £ Xs, a ,2 of the singular 
initial value problem. 



Proof of Proposition 3.5' We define the operator G :=Mo F. Here H is the solution operator as 
in Proposition 3.3 and F is the source-term operator as before. We find that under the hypothesis, 
the operator G is a contraction on closed bounded subsets of Xs^ a .2 if (5 is a sufficiently small. 
Hence the iteration sequence defined by w/j+i = G[wj] for J > 1 and, say, Wi = converges to a 
fixed point w £ Xs, a ,2 with respect to the norm || • \\£ a 2 . Because of the properties of H, a fixed 
point of G is a solution of the SIVP. Hence, we have shown existence of solutions. Uniqueness 
can be shown as follows. Given any other solution w in Xs a 2, it is a fixed point of the iteration 
Wj + \ = G[wj]. Because G is a contraction, there, however, only exists one fixed point, and hence 
w = w. □ 

For the analogous result of infinite differentiability, we only need to substitute the Lipschitz 
continuity property in the previous result as follows. For a constant e > 0, every sufficiently 
small S > and every non-negative integer fc, the operator F maps Xg , a ,k+i into Xs iU +e,k and, 
moreover, for each r > 0, there exists C > (independent of S) so that 

\\F[w] - F[w]\\ s ^ k <C\\w- w\\Z a , k+1 (3.7) 



for all w,w £ B r (0) n Xg, a ^k+i C Xs ta ,k+i- Then, there exists a unique solution w £ Xs t a,oa 
of the singular initial value problem. In order to prove this result, we first proceed as in the 
previous proposition for finitely many derivatives. The only remaining task is to show that for 
k — > co, we are allowed to choose some non- vanishing S. This can be done with a standard 
argument for hyperbolic equations. The set B r (0) is defined with respect to the norm || • \\£ a k+1 . 

We note that the constant C is allowed to depend on k and is not required to be bounded for 
k — > oo. Note that the Lipschitz estimate involves the norm || • \\J a k+1 , while the elements for 
which this estimates needs to be satisfied are required to be only in the subspace Xs, a ,k+i of 
Xs.a,k+i- The main advantage of the C°° result over the finite differentiability case is that we 
only need to check that F maps Xg, a ,k+i m to Xs^ a +n,k, instead of Xs, a ,k+i into Xs, a +e,k- 

3.3 Standard singular initial value problem 

The following discussion is devoted to the case when the function u is given by the canonical 



two-term expansion (2.6). In this case, we speak of the standard singular initial value 
problem. 
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Theorem 3.6 (Well-posedness of the standard singular initial value problem in Xs.a,2)- Given 
arbitrary asymptotic data G H 3 (U), the standard singular initial value problem admits 

a unique solution w G X$ i0lt 2 for a, 5 > 0, provided S is sufficiently small and the following 
conditions hold: 

1. Positivity condition. Suppose that we can choose a > so that the energy dissipation 
matrix (3.3 1 is positive definite at each (t,x) € (0,6) x U for a constant r\ > 0. 

2. Lipschitz continuity property. For the given a > 0, the operator F satisfies the Lipschitz 
continuity property stated in Proposition 3.5 for all asymptotic data G H 3 (U) for 



some e > 0. 

3. Integrability condition. The constants a and e satisfy 

a + e<2(f3(x) + l)-R(X 1 (x)-X 2 (x)), x G U. (3.8) 

We note that the regularity assumptions on the asymptotic data can certainly be improved. 
An analogous theorem can be formulated for the C°°-case. 



Proof. We can apply Proposition 3.5 if we are able to control the additional contribution of the 
term L[u] which has to be considered as part of the source-term. It has no contribution to the 
Lipschitz estimate (3.6), but we have to guarantee that under these hypotheses, L[u] G Xg^ a +e,i 
for the given constant e. This is indeed the case if (3.8) holds. □ 

Example 3.7. Consider the second-order hyperbolic Fuchsian equation 

D 2 v - XDv - t 2 d 2 x v = 0, 
with a constant X. This is the Euler-F 'oisson-Darboux equation. In the standard notation it is 

d 2 t v-dlv=\{\-l)d t v. 

Note that A = 1 is the standard wave equation, and in this case, the standard singular initial 
value problem reduces to the standard Cauchy problem. 

1. Case A > 0. With our notation, we have \\ = 0, A2 = —A, [3 = 0, v = 1 and f = 0. Thus 
the leading-order term for the standard singular initial value problem is 



u(t, x) 



u*(x) + u**(x)t x A > 0, 
u* (x) log t + u. JfSf (x) A = 0. 



The positivity condition of the energy dissipation matrix (3.3) is satisfied precisely for 
a > 1 — A and all sufficiently small 77 > 0. The integrability condition (3.8) is satisfied 
precisely for A < 2 — a. Hence, our previous proposition implies that the singular initial 
value problem is well-posed, provided 

< A < 2. (3.9) 

Namely, in this case there exists a solution w in Xg a 2 for some a > for arbitrary 
asymptotic data in H 3 (U). 
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2. Case A < 0. With our notation, we have Ai = |A|, A2 = 0, f3 = 0, v = 1 and f = 0. TTie 



positivity condition of the energy dissipation matrix (3.3) is satisfied precisely for a > 1 — |A| 
and aH sufficiently small 77 > . TTie integrability condition (3.8 I is satisfied precisely for 
|A| < 2 — a. Hence, our previous proposition implies that the singular initial value problem 
is well-posed, provided 

-2 < A < 0. 

Namely, in this case there exists a solution w in Xs. a ,2 for some a > for arbitrary 
asymptotic data in H 3 (U). 

It turns out that general smooth solutions to the Euler-Poisson-Darboux equation can be 
expressed explicitly by a Fourier ansatz in x and by Bessel functions in t. It is then easy to check 



that (3.9) (and similarly for A < 0) is sharp: While for < A < 2, all solutions of the equation 
behave consistently with the two-term expansion at t = 0, this is not the case for A > 2 for 
general asymptotic data. Hence the standard singular initial value problem is not well-posed for 
A > 2. This is completely consistent with our heuristic discussion in Section [2~3| underlying the 
canonical guess for the leading-order term. If e.g. A = 2, the assumption that the source-term 
t 2 d^v is negligible at t — fails since it is of the same order in t at t — as the second leading- 
order term. However, we can see in the proof of Theorem |3.6| that in the special case u* = 



(and arbitrary it**), the integrability condition (3.8) can be relaxed. For this special choice of 



data, solutions to the singular initial value problem exist even for A > 2. 

It turns out that, often in applications, the three conditions in Theorem |3.6| cannot be satisfied 
simultaneously While it is often possible to find constants a and e in accordance with the second 
and third condition, it can turn out that the corresponding choice of a is too small to make the 
energy dissipation matrix positive definite. In order to circumvent this problem, the trick is 



to choose canonical expansions of higher order Vj, see Section 2.3 as the leading-order term u 
for sufficiently large j. We refer to the singular initial value problem based on this choice of 
leading-order term as singular initial value problem of order j. For j = 1, it reduces to the 
standard singular initial value problem; hence we will focus on the case j > 2 in the following. 
Note that, if w is a solution of the singular initial value problem of order j, it is also a solution 
of the standard singular initial value problem. However, if there is only one solution w of the 
singular initial value problem of order j for given asymptotic data, it does not necessarily mean 
that w is the only solution of the standard initial value problem for the same asymptotic data. 

For the statement of the following theorem, we need the following notation. For all w € Xg a h 
(or w <E Xg, a ,k respectively), we introduce the functions Es, a ,k[w] ■ (0,(5] — > K (or Es,a,k[ w ] '■ 
(0, S] — > M respectively) which are defined in the same way as the respective norms, but the 
supremum in t has not been evaluated yet. In particular, this means that Es ia ,k[w] (or Es ia ,k[w]) 
is a bounded continuous function on (0, 5]. 

Theorem 3.8 (Well-posedness of the singular initial value problem of order j in Xg >a p). Given 
any integer j > 2 and any asymptotic data u*,u +>t € H mi (U) with m\ = 2j + I, there exists 
a unique solution w G Xg, a ,2 of the singular initial value problem of order j for some a > 0, 
provided 

1. F maps Xs,a, mi into Xs^+e.mt-i for all asymptotic data it*, it** £ H mi (U) for some e > 
and a given by a := a + (j — 2)ne. for an arbitrary k < 1. 

2. The characteristic speed satisfies 

2((3(x) + 1) > Ke for allx eU 
for the same constant k chosen earlier. 
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3. F satisfies the following Lipschitz condition: for each r > there exists a constant C > 
(independent of 5) so that 



E s ,a+e,i[F[w} - F[w}}(t) < CE Si& , 2 [w - w](t) 



for all t £ (0, S] and for all w, w £ B r (0) C Xg ^,2 



4- The energy dissipation matrix (3.3) (evaluated with a) is positive definite at each (t,x) £ 
(0, 8) x U for a constant r\ > 0. 

The third condition above is meaningful since both sides of the inequality are continuous and 
bounded functions on (0, 6]. Note that this theorem can be formulated without difficulty for the 
C°°-case and indeed leads to a simpler statement. 

In effect we have obtained a value of a which increases with j and henceforth improves the 
positivity of the energy dissipation matrix. The main prize to pay here is that the asymptotic 
data must be sufficiently regular and that we must live with a loss of regularity in space. 



3.4 The Fuchsian numerical algorithm 

We proceed with the numerical implementation of our approximation scheme. For linear source- 
terms we have shown that the solution of the singular initial value problem can be approximated 
by solutions to the regular initial value problem. We have established an explicit error estimate 
for these approximate solutions. For the nonlinear case, an additional fixed point argument was 
necessary for the proof, but the Lipschitz continuity condition should allow us to extend the 
error estimates to nonlinear source-terms. 

The regular initial value problem for second-order hyperbolic equations corresponds to the 
standard initial value problem of a system of (nonlinear) wave equations with initial time to > 0, 
and there exists a huge amount of numerical techniques for computing solutions 30 , 33! However, 
a second-order Fuchsian equation written out with the standard time-derivative dt (instead of 
D) clearly involves factors l/t or 1/t 2 . Although these are finite for the regular initial value 
problem, they still can cause severe numerical problems when the initial time t$ approaches 
zero, due to the finite representation of numbers in a computer. In order to solve this problem, 
we introduce a new time coordinate r := logt, and observe that D = d T . For instance, the 
Euler-Poisson-Darboux equation becomes 

d 2 v — Xd T v — e 2r d 2 v = 0, 

where v is the unknown and A is a constant. We have achieved that there is no singular term in 
this equation; the main price to pay, however, is that the singularity t — has been "shifted to" 
t = — oo. Another disadvantage is that the characteristic speed of this equation (defined with 
respect to the r-coordinate) is e T and hence increases exponentially with time. For any explicit 
discretization scheme, we can thus expect that the CFL-conditiorj^] is always violated from some 
time on. We must either adapt the time step to the increasing characteristic speeds in r, or, 
when we decide to work with a fixed time resolution, accept the fact that the numerical solution 
will eventually become instable. However, this is not expected to be a severe problem since 
one can compute the numerical solution with respect to the r-variable until some finite positive 
time when the numerical solution is still stable and then, if necessary, switch to a discretization 
scheme based on the original i-variable. For all the numerical solutions presented in this paper, 
however, this was not necessary. 

2 The Courant-Friedrichs-Lewy (CFL) condition for the discretization of hyperbolic equations with explicit 
schemes is discussed, for instance, in I30| . 
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We can simplify the following discussion slightly by writing (and implementing numerically) 
the equation not for the function v but for the remainder w = v — u. Henceforth we assume 
that u is the canonical two-term expansion determined by given asymptotic data. We have to 
solve the equation for w on a time interval [r , 5] for some r € K successively going to — oo with 
regular data 

w(tq,x) = 0, d T w(TQ,x) = 0, x £U. 

Inspired by Kreiss et al. in [3T] and by the general idea of the "method of lines" , see [30] , we 
proceed as follows to discretize the equation. First we consider second-order Fuchsian ordinary 
differential equations (written for a scalar equation now for simplicity) 

d 2 T w + 2ad T w + bw = /(r), 

where / is a given function and the coefficients a and b are constants. We discretize the time 
variable r so that r n := To + nAr, w n := w(T n ) and /„ := f(r n ) for some time step At > and 
n G N. Then the equation is discretized in second-order accuracy as 

(Arp + 2Ar + " = U (3 0) 

Solving this for w„+i allows to compute the solution w at the time t„ + i from the solution at the 
given and previous time t„ and r„_i, respectively. At the initial two time steps To and t\, we 
set, consistently with the initial data for w at To above, 

w a =0, Wl = ^(At) 2 /(t ). (3.11) 

We will refer to this scheme as the Fuchsian ODE solver. 

The idea of the method of lines for Fuchsian partial differential equations is to discretize also 
the spatial domain with the spatial grid spacing Ax and to use our Fuchsian ODE solver to 
integrate one step forward in time at each spatial grid point. The source-term function /, which 
might now depend on the unknown itself and its first derivatives, is then computed from the data 
on the current or the previous time levels. Here we understand that spatial derivatives are part 
of the source-term and are discretized by means of the standard second-order centered stencil 
using periodic boundary conditions. A problem is that /, besides spatial derivatives, can also 
involve time derivatives of the unknown w (in fact this can also be the case for Fuchsian ordinary 
differential equations when the source term depends on the time derivative of the unknown) . In 
order to compute those time derivatives in second-order accuracy without changing the stencil 
of the Fuchsian ODE solver, we made the following choice. In the code we store the numerical 



solution not only on two time levels, as it is necessary up to now for the scheme given by (3.10) 



and (3.11), but on a further third past time level. The time derivatives in the source-term can 
then be computed in second-order accuracy from data exclusively at the present and previous 
time steps as follows 

d T w{T n ) = — + 0((At) ). 

For this, we need to initialize three time levels at r = To and hence we set u>2 — 2(At) 2 /(tq), in 



addition to (3.11 1 



3.5 An example: the Euler-Poisson-Darboux equation 

We present numerical test results for the Euler-Poisson-Darboux equation now. Recall from 



Example 3.7 that the singular initial value problem with two-term asymptotic data for this 



1G 



equation is well-posed in particular for < A < 2, and in general becomes ill-posed for A > 2. 
The singular initial value problem for A > considers solutions of the form 

v(t, x) — u* (x) + u**(ar)t + w(t, x), 

with remainder w. For the purposes of this test, we choose the asymptotic data = cosx, 
m** = 0. Note that in this case, this leading-order behavior is consistent even with the case 
A = 0. But according to our previous discussion, it is not consistent with A = 2, and we expect 
that this becomes visible in the numerical solutions. For it** = and < A < 2, we can show 
that the leading-order behavior of the remainder at t = is 

^ «) = {-W^Yf + 8(2-AK4-A) ^ + ' • ') ' (3 - 12) 




Figure 1: Numerical solutions Euler-Poisson-Darboux equation (as explained in the text). 

First we confirm that the numerical solutions converge in second-order when At and Aa; 
are changed proportionally to each other for a given fixed choice of initial time tq > 0. In the 
following, we choose the resolution so that discretization errors are negligible relative to other 
errors. In Figure [T] now, we show the following results obtained with N — 20, At = 0.003. 
Here, N is the number of spatial grid points, i.e. one has Ax = 2w/N. The CFL-parameter 
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is At/A:r « 0.01, and we find that the runs are stable for all r < 5. For each of the plots of 
Figure [I] we fix a value of A and study the convergence of the approximate solutions to the 
(leading-order of the) exact solution (3.12) for various values of the initial time To- We plot the 
value at one spatial point x = only. The convergence rate for r — > — oo is fast if A = 1 or 
A = 0.01, but becomes lower, the more A approaches the value 2, where it becomes zero. This is in 



exact agreement with our expectations and consistent with the error estimates in Proposition 3.3 
Hence the numerical results are very promising and confirm the analytic expectations. 

Let us comment on numerical round-off errors. All numerical runs in this paper were done 
with double-precision (binary64 of IEEE 754-2008), where the real numbers are accurate for 16 
decimal digits. However, for the case To = —20 for instance, the second spatial derivative of 
the unknown in the equation is multiplied by exp(— 40) w 10 -18 at the initial time which is not 
resolved numerically and hence could possibly lead to a significant error. This, however, does 
not seem to be the case since we obtained virtually the same numerical solution with quadruple 
precision (binaryl28 of IEEE 754-2008), i.e. when the numbers in the computer are represented 
with 34 significant decimal digits. 



4 Application to Gowdy spacetimes 
4.1 Background material 

Let us provide some background material on Gowdy spacetimes [231 115j . Introduce coordinates 
(t,x, y, z) such that (x,y,z) describe spatial sections diffeomorphic to T 3 while t is a timelike 
variable. We can arrange that the Killing fields associated with the Gowdy symmetry coincide 
with the coordinate vector fields d y , d z in a global manner so that the spacetime metric reads 

g = e A/2 {-dt 2 + dx 2 ) + t {e p (dy + Qdzf + e~ p dz), t > 0. 

Hence, the metric depends on three coefficients P = P(t,x), Q = Q(t,x), and A = A(t,x). We 
also assume spatial periodicity with periodicity domain U := [0,27r). 

In the chosen gauge, the Einstein's vacuum equations imply the following second-order wave 
equations for P,Q, 



)tt H — - — Qxx — — 2(PtQt — PxQx), 



(4.1) 



which are decoupled from the wave equation satisfied by the third coefficient A: 

A« - Kxx = Pi ~ i? + e 2P (Q 2 x - Q 2 ). (4.2) 

Moreover, the Einstein equations also imply constraint equations, which read 

A, = 2t (P x P t + e 2P Q x Q t ) , (4.3a) 

A t - * (Pi + te 2P Q 2 x + Pi + e 2P Q 2 ) . (4.3b) 



It turns out that (4.2 1 can sometimes be ignored in the following sense. Given a time to > 0, 
we can prescribe initial data (P,Q)\t f° r the system (4.1) while assuming the condition 

r-2n 



\ (P x P t + e 2P Q x Q t ) dx = at t = t . 
Jo 
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Then, the first constraint (4.3a) determines the function A at the initial time, up to a constant 



which we henceforth fix. Next, one easily checks that the solution (P,Q) of (4.1 1 corresponding 



to these initial data does satisfy the compatibility condition associated with (4.3) and, hence, 
(4.3) determines A uniquely for all times of the evolution. Moreover, one checks that (4.2) is 



satisfied identically by the constructed solution (P,Q,A). One can also consider the alternative 
viewpoint which follows from the natural 3 + 1-splitting and treats the three equations (4.1 )-(4.2 ) 



as an evolution system for the unknowns (P, Q,A), and (4.3) as constraints that propagate if 
they hold on an initial hypersurface. In any case, equations (4.1) represent the essential set 
of Einstein's field equations for Gowdy spacetimes. We refer to (4.1) as the Gowdy equations 
and focus our attention on them in most of what follows. An alternative, more geometrical 
formulation of Einstein's field equation for Gowdy symmetry has been introduced in [3] . 



4.2 Heuristics about singular solutions of the Gowdy equations 

We provide here a formal discussion which motivates the (rigorous) analysis in subsequent sec- 
tions. Based on extensive numerical experiments [6, 5, 3], it was first conjectured (and later 
established rigorously [39l [40]) that as one approaches the singularity the spatial derivative of 
solutions (P, Q) to (4.1 ) becomes negligible and (P, Q) should approach a solution of the ordinary 
differential equations 



Qt 



-2P t Qt 



(4.4) 



These equations are referred to, in the literature]^] as the velocity term dominated (VTD) equa- 
tions. Interestingly enough, they admit solutions given explicitly by 



P(t,x) = log (at k (l + C 2 t- 2k )), Q(t,x) = £ - 



-2k 



a{l + Ct- 2k ) 



(4.5) 



where x plays simply the role of a parameter and a > 0, £, £, k are arbitrary 2-7r-periodic functions 
of x. 



Based on (4.5), it is a simple matter to determine the first terms in the expansion of the 



function P near t = 0, that is 



lim V = limtP t (t,x) = -\k\, 
i-s-0 lost t->o y ' ; 11 



hence 



lim (P(t, x) 
t->o v 



log a, 

\k(x)\\ogt) = <p(x), cp := { l g(a(l + C 2 )), 

log« 2 ), 



k < 0, 
k = 0, 
k > 0. 



3 To our knowledge, this terminology has been introduced in |20| and, in the context of Gowdy spacetimes, was 
first used in [27]. 
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Similarly, for the function Q we obtain if £ ^ 0, 



lim Q(t, x) = q(x), 



c 



fc < o, 

-, k = 0, 



c 

a 



fc < 0, 

lim tr 2 ^ (Q(t, x) - q{x)) = tf>(x), ip := { 0, fc = 0, 
t— >o ] 



If C = 0, then g = C- 



From (4.5), we thus have the expansion 
P = -\k\ log* + 95 + 0(1), 



Q^ + t^V + o^ 1 * 1 ), 



2|A+ 



(4.6) 



in which fc, cp, q, %p are functions of x. In general, P blows-up to +00 when one approaches 
the singularity, while Q remains bounded. Observe that the sign of k is irrelevant as far the 
asymptotic expansion is concerned, and we are allowed to restrict attention to fc > 0. 



By plugging the explicit solution into the nonlinear terms arising in (4.4) one sees that e 2P Q 



is of order i 2 d fc l _1 ) which is negligible since the left-hand side of the P-equation is of order i~ 2 , 
at least when fc =^ 0. On the other hand, the nonlinear term PtQt is of order t 2 ^' 1 ', which is 
the same order as the left-hand side of the Q-equation. It is not negligible, but we observe that 
PtQt has the same behavior as —(\k\/t)Q t . 



In fact, observe that the homogeneous system deduced from (4.4) 



Pti 



P 



= 0, 



1 — 2k 

Qu + — — Qt = 0, 



(4.7) 



is solved precisely by the leading-order terms in (4.6). This tells us that, as t — > 0, the term 
e 2P Q 2 is negligible in the first equation in (4.4), while PtQt + (\k\/t)Q t is negligible at t = 0. This 



discussion hence allows us to conclude that as far as the behavior at the coordinate singularity 



t = is concerned, the nonlinear VTD equations (4.4) are well approximated by the system 



(4.7) 



We return now to the nonlinear terms which were not included in the VTD equations, but 



yet are present in the full model (4.1). Allowing ourselves to differentiate the expansion (4.6), 
we get the following leading-order terms at t = 0: 

' t -2\k\ e 2 Vq 2 + _ ^ 

e 2P Q 2 x = { 2e 2 ^\k\ x ^\ogt+. 



Px Qx — 



3 2 ^ + ..., 

-log* \k\ x q x + 

(f x Qx + 

-21og 2 « 2 l fc l(|fcU)V 



q x = 0, 
q x = 0, 



\k\ x ¥=o, 
\k\ x = 0, 



\k\ x ,q x ^ 0, 
|fc| x =0, (p x ,q x ^0, 

\k\ x = q x = 0, (f x ,ijj x ^0. 



To check (formally) the validity of the expansion (4.6) we now return to the full system 



Consider the nonlinear term e 2P Q 2 in (4.1), and observe the following 
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Case q x 7^ everywhere on an open subinterval of [0, 2it\. Then, on the one hand, the 
left-hand side of the first equation in (4.1 ) is of order t~ 2 , at most. On the other hand, the 



term e 2P Q 2 is negligible with respect to t 2 if and only if \k\ < 1 and is of the same order 
if = 1. 

• Case q x = on an open subinterval of [0, 2ir]. Then, e 2P Q 2 X is negligible with respect to 
t~~ 2 , and no condition on the velocity k is required on that interval. 

• Case q x [xq) — at some isolated point Xq. Then, no definite conclusion can be obtained 
and a "competition" between \k\ (which may approach the interval [0, 1]) and q x (which 
approaches zero) is expected. 

Similarly, at least when \k\ x q x ^ 0, the nonlinear term P X Q X is of order \ogt and, therefore, 
negligible with respect to t 2 d fe l -1 ) (given by the left-hand side of the second equation in (4.1 )) if 



and only if \k\ < 1. Points where \k\ x or q x vanish lead to a less singular behavior and condition 
on the velocity can also be relaxed. 

The formal derivation above strongly suggests that we seek solutions to the full nonlinear 



equations admitting an asymptotic expansion of the form (4.6 1, that is 

P = -k log* + ip + o(l), Q = q + t 2k {ijj + o(l)), 

where k > and ip, q, %j> are prescribed. In other words, these solutions asymptotically approach 
a solution of the VTD equations and, in consequence, such solutions will be referred to as 
asymptotically velocity term dominated (AVTD) solutions [27] , 

Based on this analysis and extensive numerical experiments, it has been conjectured that 
asymptotically as one approaches the coordinate singularity t — the function P(t, x) / \ogt 
should approach some limit k — k(x), referred to as the asymptotic velocity, and that k(x) 
should belong to [0, 1) with the exception of a zero measure set of "exceptional values". The 
reason for this name of k is the following. Based on the work by Geroch [25] , it was noted 



by Moncrief [32] that the evolution equations (4.1) for P and Q can be considered as wave map 



equations with the hyperbolic space as the target space. If a solution of these equations has an 



expansion of the form (4.6) at t = 0, then the velocity of the image points of this map, which 
must be defined with the correct convention of the sign and must be measured with respect to 
the hyperbolic metric, approaches k as t — > 0. 

It was demonstrated in [5] that solutions of the Gowdy equations which are compatible with 



(4.6) approach certain Kasncr solutions at t = 0, with possibly different parameters along each 



timeline to the singularity. 

4.3 Gowdy equations as a second-order hyperbolic Fuchsian system 



The first step in our (rigorous) analysis of the Gowdy equations ( |4.1| now is to write them as a 
system of second-order hyperbolic Fuchsian equations. After multiplication by t 2 , the equations 



(4.1) immediately take the second-order hyperbolic Fuchsian form 

D 2 P = t 2 dlP + e 2P (DQ) 2 - t 2 e 2P {d x Qf 
D 2 Q = t 2 d 2 Q - 2DPDQ + 2t 2 d x Pd x Q. 

The general canonical two-term expansion then reads 

P(t, x) = F* (x) log t + P**(x) + ... 
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for the function P and, similarly, an expansion Q*(x) logt + Qt,*(x) + . . . for the function Q with 
prescribed data Q* , Q** . At this stage, we do not make precise statements about the (higher- 
order) remainders, yet. In any case, the Fuchsian theory does not apply to this system due to 
the presence of the term —2DPDQ (with the exception of the cases P» = or =0). Namely, 
this term does not behave as a positive power of t at t = when we substitute P and Q by their 
canonical two-term expansions, but this is required by the theory. The reason for this problem 



is the significance of the nonlinear term in the source-term as found in Section 4.2 cf. (4.7). 

We propose to add the term ~2kDQ to the equation for Q where A; is a prescribed (smooth, 
spatially periodic) function depending on x, only. The function k will play the role of the 
asymptotic velocity mentioned before. This yields the system of equations 



D 2 P = t 2 d 2 x P + e 2P (DQ) 2 - t 2 e 2P (d x Q) 2 , 
D 2 Q - 2kDQ = t 2 d' 2 Q - 2(fc + DP)DQ + 2t 2 d x Pd x Q. 



(4.8) 



The resulting system is of second-order hyperbolic Fuchsian form with two equations, corre- 
sponding to 



A« = A« = 0, 



a! 2) = 0, 



A^ ^ — 2k, 



Here, the superscript determines the respective equation of the system (4.8). If we assume that 



A: is a strictly positive function, as we will do in all of what follows, the expected leading-order 
behavior at t = given by the canonical two-term expansions is 



P(t, x) = P* (x) log t + P« (x) + 
Q{t,x) = Q*(x) + Q**(x)t 2k ( x U 



(4.9) 



One checks easily that the problem associated with the term —2DPDQ before does not arise if 
P* = —k. Indeed, the canonical two-term expansion (4.9) is consistent with the heuristics of the 



Gowdy equations above and we recover the singular initial value problem studied rigorously in 
[28l 136] and numerically in pQ . We only mention here without further notice that the case k = 
with the logarithmic canonical two-term expansion for Q is covered by the following discussion. 
Furthermore, the case of k vanishing at only certain points may be also included via a suitable 



renormalization of the asymptotic data, see (2.7) 
When P, = 



k, the function k plays a two- fold role in (4.8). On the one hand, it is an 



asymptotic data for the function P and, on the other hand, it is a coefficient of the principal 
part of the second equation. In order to keep these two roles of k separated in a first stage, we 
consider the system 



D 2 P = t 2 d 2 x P + e 2P (DQ) 2 - t 2 e 2P (d x Q) 2 , 
D 2 Q - 2kDQ = t 2 d x Q - 2(-P* + DP)DQ + 2t 2 d x Pd x Q, 



(4.10) 



instead of (4.8). Studying the singular initial value problem with two-term asymptotic data 



means that we search for solutions to (4.10) of the form (as t 0) 



P(t, x) = P* (x) log t + P** (x) + w (1) (t, x) , 
Q(t,x) = Q m {x) + Q**(x)t 2k ^ +w {2) (t,x), 



(4.11) 



,,, and remainders w^K After studying the 



for general asymptotic data P*, P**, Q*, Q 
well-posedness for this problem, we can always choose P* to coincide with — k and, therefore, 



recover our original Gowdy problem (4.8)-(4.9). For simplicity in the presentation, we always 
assume that k is a C°° function. 
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In the following discussion, we write the vector-valued remainder as w := (w/ 1 ', w^), and we 
fix some asymptotic data P*, P**, Q*, and Q** and choose the leading-order term u according 
to (4.11). The source-term operator F[w](t,x) =: (Fi[w](t, x) , F 2 [w](t, x)) reads 

Dw^y 



Fi[w] 



(2kt 2k Q 



(t p - e p " e w<1) (td x Q* + 2d x kt 2k t logiQ** + t 2k td x Q„ +td x w 



(2)- 



and 
F 2 [w] 

= - 2Dw (1) (2kt 2k Q„ + Dw {2) ) 

+ 2(td x P*\ogt + td x P„+td x w {1) ^(td x Q*+2d x k^ 

4.4 Properties of the source-term operator 

To establish the well-posedness of the singular initial value problem for the Gowdy equations, 
we need first to derive certain decay properties of the source-term operator F consistent with 



Section 3 Let us introduce some notation specific to the Gowdy equations. Let X 



be the 

space denned as above based on the coefficients of the first equation in (14. 101) and, similarly, 



let X 



(2) 



(1) 

(5,Qi ,fc 



be the space associated with the second equation. By definition, a vector-valued 



map w — (toM, to( 2 )) belongs to Xs, a ,k precisely if G X s ^ fe and G Xg ^ fe , with 
o := (ax, a 2 ). An analogous notation is used for the spaces Xg 1 ^ k , Xjj^ k and Xs, a ,k- 



Now we are ready to state a first result about the source-term of (4.10) 



Lemma 4.1 (Operator F in the finite differentiability class). Fix any 5 > and any asymptotic 
data P* , P** , Q* , <3** G H m (U), m > 2. Suppose there exist e > and a continuous function 
a = (0:1,012) : U — > (0,oo) 2 so that, at each x G U, 



ai(x) + e < min (2(P*(x) + 2k(x)), 2(P*(x) + 1)), 
a 2 (x) +e < 2(1 - k(x)), 
a\(x) — a 2 (x) > e + min (0,2fe(x) — l), 
e < 1. 



(4.12a) 
(4.12b) 
(4.12c) 
(4.12d) 



Then, the operator F associated with the system (4.10) and the given asymptotic data maps 
Xs, a ,m into Xs^ a +t,m-i o,nd satisfies the following Lipschitz continuity condition: For each r > 
and for some constant C > ( independent of 8 ), 



E 5 , 



for all w,w G B r C X$ 



F[w]-F[w] (t)<CE s , a , m [w-w](t), te(0,S\ 



where B r denotes the closed ball centered at the origin. 



In this lemma, since P* G if 1 (17), in particular, a standard Sobolev inequality implies that P* 
can be identified with a unique bounded continuous periodic function on U, and the inequality 



(4.12a) makes sense pointwisc. 



Proof. Consider the expression of F given before. Let w G Xs. a , m for some (so far unspecified) 
positive spatially dependent functions 011,02, hence G X/ja m an d G Xg^ m . By a 
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standard Sobolev inequality (since m > 2 and the spatial dimension is 1), we get that F[w\(t, ■) £ 
H m ~ 1 (U) for all t £ (0,8). Namely, if m > 2 we can control the nonlinear terms of F[w](t,-) in 
all generality for a given t > if any factor in any term of F[w](t, •), after applying up to m — 1 
spatial derivatives, is an element in L°°(U) - with the exception of the mth spatial derivative of 
w which is only required to be in L 2 (U). This is guaranteed by the Sobolev inequalities. Having 
found that F[w]{t, •) £ H 171 - 1 ^) for all t £ (0,<5], it is easy to check that Fi[w] £ x£l i+e if 



ax(x) +e < min (2(P*(x) + 2k(x)),2{P*(x) + 1)). 



x£U. 



(4.13) 



Even more, condition (4.13) implies that D l Fi[w] £ jj 1 ] ^ for all I < m — 1 



Considering now spatial derivatives, we have to deal with two difficulties. The first one is 



that logarithmic terms arise with each spatial derivative. We find d x D l Fi[w] £ xj- 1 ^ , n 
I < m — 1 and k < m — 2 and k + I < m — 1 (excluding first the case k = m — 1, I 



for all 



Oil 



(x) + e < min (x) + 2k(x)), 2(P*(x) + 1) 



.t € 17. 



0) provided 
(4.14) 



A second difficulty arises in the case k = m — 1, I = 0. Namely, since w £ Xs, a ,m (and 
not in X s>a<rn ), it follows that in particular td™wW ~ i 2fc+Q 2 ( anc i not t i+2fe+a2). note ^hat the 
function /3 which determines the behavior of the characteristic speeds at t = is identically zero 
in the case of the Gowdy equations. The potentially problematic term is hence of the form AB 
with 



A :=t 
B:=t p *e p - 



(td x Q* + 2d x kt 2k t log tQ„ + t ZK td x Q^ + td x w {2 >) 



2k a 



>e wil) (d^itdzQ* + 2d x kt 2k t logtQ^ + t 2k td x Q„) + td^w^) 



originating from taking m — 1 spatial derivatives of F\ [w] . To ensure d™ 1 Fi[w] ^ Xg 1 ^ i+e Q , we 
need 

ai (x) + e< (P,(x) + l) + (P,(x) + 2k(x)+a 2 (x)), x £ U. (4.15) 



If (4.141 is satisfied, we have (for all x) 

L (x) + e < min (2(P*(x) + 2k(x)), 2(P*(as) + 



«i( 



< (P m (x) + l) + (P m (x)+2k(x)) 



and, thus, (4.151 follows from (4.14|. In conclusion, (4.141 is sufficient to guarantee that Fi[w] £ 



X 



(i) 

<5,ai+e,m — 1 ' 

Let us proceed next with the analysis of the term F 2 [w] . If 



ai(x) - a 2 (x) > e, a 2 (x) + e < 2(1 - k(x)), x £ U, (4.16) 

(2) (2) 

then [iw] £ Xg c , 2+c - This inequality also implies that all time derivatives are in Xg ^ 2+e as 
before. We have to deal with the same two difficulties as before when we consider spatial deriva- 
tives of P2[f]- O n the one hand, equality in (4.16) cannot occur due to additional logarithmic 
terms. On the other hand, we must be careful with the (m — l)-th spatial derivative of F 2 [w]. 
Here, the two problematic terms are of the form AB with either 



A :=d?- 1 (td x P, log* + td x P, 
B :=td x Q* + 2d x kt 2k t\ogtQ„ 



t tdxQ * 



td x 



,(2) 
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or else 



A :=td x P*\ogt 
B:=d™-\td x Q 



td x P,,+td x w (1 \ 
+ 2d x kt 2k t\ogtQ if 



The first one is under control provided ai(x) + 1 > 2k(x) + a2(x) + e, for all x G U, while for the 
second one it is sufficient to require e < 1. The claimed Lipschitz continuity condition follows 
from the above arguments. □ 



Positive functions ct\ and 02 and constants e > satisfying the hypothesis of Lemma 4.1 can 



obviously exist only if k(x) < 1 for all x G U (due to (4.12b)). In Lemma 4.3 below we identify a 
special case where this limitation is avoided. Hence, we make the assumption that < k(x) < 1 
for all x, which is consistent with our formal analysis in Section |4.2| As a consistency check for 
the case of interest P» = — k, let us determine under which conditions the inequalities (4.12) 



can be hoped to be satisfied at all. For this, consider (4.12a) and (4.12c) in the borderline case 
c*2 = e = 0. This leads to the condition < k < 3/4, which shows that Lemma 4.1 does not 
apply within the full interval < k < 1. It is interesting to note that Rendall was led to the same 
restriction in |36j . but its origin stayed unclear in his approach. Here, we find that this is caused 



by the presence of the condition (4.12c) in particular which reflects the fact that w is an element 
of the space Xs, a , m rather than of the smaller space X^ a>m . Interestingly, we can eliminate this 
condition and, hence, retain the full interval < k < 1, when we consider the C°°-case, instead 
of finite differentiability. See also [3] for a detailed discussion of the different intervals of k. 

Lemma 4.2 (Operator F in the C°° class. General theory). Fix any 6 > and any asymptotic 
data P., P**, Q*, Q** G C°°(U). Suppose there exist a constant e > and a continuous functions 
a = (0:1,0:2) : U — > (0,oo) 2 such that, at each x G U, 



ax(x) + e < mm(2(P*(a;) +2k(x)),2(P*(x) 
a 2 (x)+e< 2(1 - k(x)), 
cei(x) — 0:2(0;) > e. 



1)). 



(4.17a) 
(4.17b) 
(4.17c) 



Then, for each integer m > 1, the operator F maps Xs^ a ,m into Xs ta + t ^ m -i and satisfies the 
following Lipschitz continuity property: for each r > and some constant C > (independent 
ofS), 

£ jia+e , m -i[FH-f[ffl]](t) <CE 5 , a<m [w-w\{t), te(0,6], 
for all w,w G B r n X StCX+ ^ m C X^ a+em . 



The proof is completely analogous to that of Lemma 4.1 Since only spaces Xs. a ,k need to 



be checked (i.e. without the tilde) in the C°°-case we obtain stronger control than in the finite 
differentiability case. In consequence, the C°°-case does not require the condition (4.12c). Thus 



k can have values in the whole interval (0, 1) as we show in detail later. In a special case, which 
will be of interest for the later discussion, however, we can relax the constraints for k even in the 
finite differentiability case. 

Lemma 4.3 (Operator F in the finite differentiability class. A special case). Fix any 8 > and 

any asymptotic data P*,P**,Q*« £ H m {U), = const, m > 2. Suppose there exist e > and 
a continuous function a = (01,02) : U — > (0,oo) 2 such that, at each x G U, 



a 1 (x) + e<2{P*{x) + 2k(x)), 

02(2) + e < 2, ai(x) — a^{x) > e 



1. 



e < 1. 



Then, the operator F satisfies the conclusions of Lemma 4-1 
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In the special case = const, we have hence characterized the map F for k being any 
positive function in the finite differentiability case. The analogous result for the C°°-case can 
also be derived. 



4.5 Well-posedness theory 



Relying on Theorem |3.6| and the results in the previous sections, we now determine conditions 
that ensure that the singular initial value problem for the Gowdy equations is well-posed. Besides 
the properties of the source-operator F already discussed, we have to check the positivity of the 
energy dissipation matrix. This leads us to the matrix 

ai -r)/2 
jV« : = | - v /2 at 





for the first component and to the matrix 

'2k 



Q.2 

-77/2 




-r?/2 

a.2 

~2d x k(tlogt) 





-2d x k{t log t) 
2k + a 2 -l 



for the second component. For the matrix to be positive, it is necessary that ai(x) > 1 for 
all x e U. However, if P* = — k, then condition (4.12a) in Lemma 4.1 in the finite differentiability 



case (or the corresponding one in Lemma 4.2 in the C°°-case) implies that ai(x) < 1. Hence, in 
the same way as in Rendall |36j . one does not arrive at a well-posedness result for the singular 
initial value problem yet. However, since the positivity of the energy dissipation matrix is the 



only part of the hypothesis in Theorem 3.6 which is is violated, we can use instead Theorem |3.8| to 
prove well-posedness of the singular initial value problem with asymptotic solutions of sufficiently 
high order j. 

Let us be specific about what we mean by j being "sufficiently large" , and we now make 
some choice for the parameters ai, a 2 and e, consistent with Lemma |4.2[ which will allow us to 
estimate the required size of j. We make no particular effort to choose these quantities optimally, 
but still the goal is to choose j "reasonably" small. Henceforth, we restrict to the C°°-case and 
P*{x) = —k(x) with < k(x) < 1 for all x G U. We introduce positive constants Mi and u 2 
(with further restrictions later) and the function x( x ) 
states that we must choose a±(x) and e so that a\{x) 



= l-2\x- 1/2] . The condition ( |4.17a[ ) 
e < x(k(x)). We set 



ai (x) := 1 - yji(k(x) - 1/2) 2 + /if, (4.18) 
and find x(fc(x)) — oti(x) > \J\ + n\ — 1 for all x g U, provided < k(x) < 1. Similarly, we set 

a 2 (x) := 1 - ^A(k(x)-l/2f + nl, (4.19) 



and it follows that a\{x) — ct2{x) > \Jl + fi 2 — \A + Mi f° r a 2 > Ml- F° r the conditions (4.17a) 
and (4.17c) to hold true, we have to choose 



< /ii < //2, and < e < min 




l,i/l 



M2 




Condition (4.17b) is then satisfied automatically. 
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Now, assume in what follows that k(x) € (1/2 — Ak, 1/2 + Ak) for all x € U for a constant 
Afc e (0, 1/2). Then it is clear that both functions ai and 02 are positive for all such k(x) if and 
only if 

Mi < M2 < v/l-4(Afc) 2 . 



This assumption will be made in the following. In Theorem |3.8[ we could choose j as small as 
possible if we pick the maximal allowed value for e. Hence, we set 




i4-i,ji + i4- ji 




We find easily that 
provided 



1 



Mi < j(mI + 2^1 + -2), 



and check that this is consistent with the condition < Mi < f-2 made before. In order to make 
a specific choice, we assume this inequality for Mi and hence obtain that 



l + /i?-l. (4.20) 
Now, in order to make the energy dissipation matrix positive, we must choose j so that for all 

xeu, 

5i(x) := cki(x) + (j — 2)ne > 1, 

ot2{x) :— a^{x) + (j — 2)Ke > 1 — 2k(x); 

cf. Theorem |3.8| These two inequalities are satisfied for all functions k under our assumptions if 
in particular 



In any case, we choose the maximal value for mi 



Ml 



^■=\\ll>i (4.221 



since this minimizes the value on the right side of (4.21). We find that for this value of mi, the 



right side of (4.211 is monotonically decreasing in M2 and diverges to +oo for M2 — > for all 
values of Afc. 

Theorem 4.4 (Well-posedness theory for the Gowdy equations). Consider some asymptotic 
data = —k, P**, Q*, Q** S C°°(U), where k is a smooth function U — > (1/2 — Afc, 1/2 + Afc) 
for a constant Ak € (0, 1/2). Then, the singular initial value problem with asymptotic solutions 
of order j has a unique solution with remainder w G ^<5.Q+(j-2)Ke,oo f or some sufficiently small 
5 > and some k < 1. Here, the exponents a — (0:1,0:2) & n d 6 ore given in ( |4.18[ ), (4.19), and 



(4.20) explicitly in terms of the data and parameters \x\ ,M2 chosen such that mi is an explicit 



expression in M2 given in (4.22) while M2 is a sufficiently close to (but smaller than) y/T— 4(Afc) 2 , 



and the order of differentiation j satisfies 

j > 2 + 



3 - 4(Afc) 2 + 2^2-4(Afc) 2 -2 
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The above condition implies that to reach Ak -> we need j > 7, while Ak — > 1/2 requires 
j ^ oo. Although our estimates may not be quite optimal, the latter implication cannot be 
avoided. 



4.6 Fuchsian analysis for the function A 



So far we have considered the equations (4.1| for P and Q. We can henceforth assume that 



these equations are solved identically for all f > (and t < 5 for some d > 0) and that hence P 



and Q are given functions with leading-order behavior (4.9 1 and remainders in a given Xs^ a ,k 



The equations which remain to be solved in order to obtain a solution of the full Einstein's 
field equations are (4.2) and (4.3). In particular we are interested in the function A in order 



to obtain the full geometrical information. We compute A from a singular initial value problem 
with "data" on the singularity analogously to P and Q. The following discussion resembles the 
previous one and we only discuss new aspects now. 



Clearly, the three remaining equations (4.2) and (4.3) for A are overdetermined, and hence 



solutions will exist only under certain conditions. Let us define the following "constraint quan- 



tities" from (4.3) 



Ci(t, x) := -d t A + t(P x ) 2 + e 2P t(Q x ) 2 + t(d t P) 2 + e 2P t(d t Q) 2 , 
C 2 (t, x) := -A, + 2P X DP + 2e 2P Q x DQ. 



Moreover, we define 



H(t,x) := -A t 



•A* 



p--Pt+e 2P {Qii-Qi) 



from (4.2). From the evolution equations for P and Q, we find the subsidiary system 



d t Ct = d x C 2 + H, 



d t C 2 = d x d. 



(4.23) 



These equations have the following consequences. Suppose that we use (4.3b) as an evolution 
equation for A. This implies that C\ = for all t > 0. Moreover, suppose that we prescribe 
data at some to > (indeed to is allowed to be zero later) so that C2{to,x) — for all x € U. 
Then the equations imply that H = and C 2 = for all t > and thus we have constructed a 



solution of the full set of field equations. Alternatively, let us use (4.2) as the evolution equation 
for A, i.e. H = 0. Suppose that we prescribe data so that Cx(to,x) = C 2 (to,x) — at some to- 
It follows that C\ = C 2 = for all t > because the evolution system (4.23) for C\ and C 2 is 



symmetric hyperbolic. Again, Einstein's field equations are solved. 

Now, we want to consider the case to = 0. First note that (4.23) is regular even at t = 0. 



Suppose that P and Q are functions with leading-order behavior (4.9 1 and remainders in a given 
Xg.a.k with k > 1. If there exists a function W3 so that 

A(t,x) = A*(a;)logt + A*,(a:) +w 3 (t,x) (4.24) 

with W3 converging to zero in a suitable norm at t = and 

A, (a) = k 2 (x), A,,(x)=Ao + 2 [ k(x)(-d x P„(x) + 2e 2P "^Q„(x)d x Q,{x)) dx, (4.25) 



where Aq is an arbitrary real constant, then, in particular, 



lim C 2 = 0. 
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It also follows that lim t _>oiCi = 0. Let us first use (4.3b I as a singular evolution equation for 



A. Since this is "only" a singular ODE, one can show easily that there exists a unique solution 
for A for t > which obeys the two-term expansion above, and hence C\ = 0. Our discussion 
before implies that H, C2 = 0. Hence we obtain a solution of the full Einstein's field equation 



for all t > 0. Alternative, choose (4.2) as the evolution equation for A now. This equation can 



be written in second-order hyperbolic Fuchsian form 



D 2 A - t 2 d 2 x A = {td x Py + (DA - (DP) 2 ) + e 2r ((td x QY - (DQ) 2 ) 



Indeed this equation is compatible with the leading-order expansion (4.24) at t — and we can 



show well-posed of this singular initial value problem in the same way as we did for the functions 
P and Q before (also going to sufficiently high-order in j). In particular, for any asymptotic 
data A* and A**, not necessarily those given by (4.25), there exists a unique solution of this 
equation A with remainder W3 in a certain space Xs t<x ,k- By means of uniqueness we find that 



the solution A of this equation coincides with the solution for A obtained using (4.3b) as the 
evolution equation. Hence, we must have H, C\, C2 = for all t > 0, and thus also this method 
yields a solution of the full Einstein's field equations. 



Note that periodicity and (4.25) implies that the asymptotic data for P and Q must satisfy 
the relation 



2,T 



k(x)(-ds:P^(x) + 2e 2P " (£) Q«(i)<9 £ Q st (£)) dx = 



for smooth solutions. 



5 Numerical experiments 

5.1 Test 1. Homogeneous pseudo-polarized solutions 

We continue our discussion with the singular initial value problem for the Gowdy equations. In 
all of what follows we consider the singular initial value problem with two-term asymptotic data 
for the Gowdy equations. As we show this works very well and we get good convergence. This 
is a strong indication that the standard singular initial value problem is well-posed. In contrast, 
recall from Theorem |4 . 4| that our analytical techniques are only sufficient to show that the initial 
value problem with asymptotic solutions of sufficiently high order is well-posed for the Gowdy 
equations. 

Before we proceed with "interesting" solutions of the Gowdy equations, let us start with a 
case for which we can construct an explicit solution and hence test the numerical implementation. 
Let P and Q be solutions of the polarized equations in the homogeneous case, i.e. set Q = and 
P(t, x) = P(t). In this case, it follows directly that the exact solution of the Gowdy equations is 

P(t) = -fclogt + P„, 

where both k and P** are arbitrary constants. The corresponding full solutions of Einstein's 
equations are Kasner solutions whose parameter^] are determined by k exclusively (P** is just a 
gauge quantity) . By a reparametrization of the Killing orbits of the form 

^2 = X2/V2 + x 3 / V2, x 3 = —X2/V2 + x 3 /v2, 

4 In the conventions of [42], we have pi = (k 2 - l)/(fc 2 + 3), p 2 = 2(1 - fc)/(fc 2 + 3), p 3 = 2(1 + fc)/(fc 2 + 3), 
and the three flat cases are realized by k = 1, fe = —1 and |fc| — > 00. 
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(a) k = 0.5. 



(b) k = 0.9. 



Figure 2: Convergence of numerical solutions of the test case 1 (as explained in the text). 



where x 2 , and x 3 are the coordinates used to represent the orbits of the polarized solution above, 
the same solution gets re-expressed in terms of functions 



P = logcosh(-fclogi + FU), Q = tanh(-fclog< + 



(5.1) 



These functions (P, Q) are again solutions of (4.1 1. Asymptotically at t — 0, they satisfy 
P=-k log t + {P„- log 2) + . . 



Q = 1 - 2e 



from which we can read off the corresponding asymptotic data. 

Now we compute the solutions corresponding to these asymptotic data numerically and com- 
pare them to the exact solution (5.1 1. We pick P** = 1, so that P„ 
— 2e~ 2 . Since the solution is spatially homogeneous 



1 - log 2, Q* = 1 and 

Q** = — 2e~ 2 . Since the solution is spatially homogeneous - in fact this is an ODE problem 
we only need to do the comparison at one spatial point. 

The results are presented in Figure [2] where we plot the difference of the numerical and the 
exact value of Q versus time for various values of To- In the first plot, this is done for k = 0.5 and 
in the second plot for k = 0.9. The plots confirm nice convergence of the approximate solutions 
to the exact solution. The fact that each approximate solution diverges from the exact solution 
almost exponentially in time is a feature of the approximate solutions themselves and not of the 
numerical discretization, as is checked by comparing two different values of Ar in these plots. 
From our experience with the Euler-Poisson-Darboux equation, we could have expected that the 
convergence rate is lower in the case k = 0.9 than in the case k = 0.5 (note that k plays the same 
role A/2). In the case of the Euler-Poisson-Darboux equation, the rate of convergence decreases 
when A approaches 2, due to the influence of the second-spatial derivative term in the equation. 
In the spatially homogeneous case here, however, this term is zero and hence this phenomenon 
is not present. The "spikes" in Figure [2] are just a consequence of the logarithmic scale of the 
horizontal axes and the fact that the numerical and exact solutions equal for some instances of 
time. 
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Figure 3: Convergence of numerical solutions of the test case 2 as explained in the text. 

5.2 Test 2. General Gowdy equations 

Now we want to study the convergence for a "generic" inhomogeneous Gowdy case (still ignoring 
the equation for the quantity A). Here we choose the following asymptotic data 

k(x) = 1/2 + Acos{x), = 1.0 + sin(x), 

P** = 1 -log2 + cos(x), Q„=-2e~ 2 , 

with a constant A e (—1/2, 1/2). We do not know of an explicit solution in this case. In Figure[3j 
we show the following numerical results for A = 0.2 and A = 0.4, respectively. For the given value 
of A, we compute five approximate solutions with initial times tq = —30, —35, —40, —45, —50 
numerically, each with the same resolution At = 0.01 and N = 80. The resolution parameters 
have been chosen so that the numerical discretization errors are negligible in the plots of Figure [3] 
Then, for each time step for r > —30, we compute the supremum norm in space of the difference 
of the remainders w^ 1 ' of the two approximate solutions given by t = —30 and r = — 35. In 
this way we obtain the first curve in each of the plots of Figure [3] The same is done for the 
difference between the cases tq = —35 and tq = —40 for all t > —35 to obtain the second 
curve etc. Hence these curves yield a measure of the convergence rate of the approximation 
scheme (without referring to the exact solution). In agreement with our observation for the 
Euler-Poisson-Darboux equation, the convergence rate is high if k is close to 1/2 and becomes 
lower, the more k touches the "extreme" values k = and k = 1. 

Much in the same way as for the Euler-Poisson-Darboux equation we find that double preci- 
sion is sufficient for these computations despite of the fact that exp(2r) is 10~ 44 for r = —50. 

5.3 Test 3. Gowdy spacetimes containing a Cauchy horizon 

The papers [THl [T71 [TH1 [Ml HZ] were devoted to the construction and characterization of Gowdy 
solutions with Cauchy horizons in order to prove the strong cosmic censorship conjecture in this 
class of spacetimes. Spacetimes with Cauchy horizons are expected to have saddle and physically 
"undesired" properties, in particular they often allow various inequivalent smooth extensions. 
This has the undesired consequence that the Cauchy problem of Einstein's field equations does 
not select one of them uniquely. Some explicit examples are known, but most of the analysis is 
on the level of existence proofs and asymptotic expansions. 
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Hence, it is of interest to construct such solutions numerically and analyze them in much 
greater detail than possible with purely analytic methods. Constructing these solutions numer- 
ically, however, is delicate since the strong cosmic censorship conjecture suggests that they are 
instable under generic perturbations. It can hence often be expected that numerical errors would 
most likely "destroy the Cauchy horizon". This is so, in particular, when the singular time at 
t = is approached backwards in time from some regular Cauchy surface at t > 0, i.e. for the 
"backward approach" . 

In the Gowdy case, where the strong cosmic censorship conjecture has been proven |40j . 
however, there are clear criteria for the asymptotic data so that the corresponding solution of 
the singular initial value problem has a Cauchy horizon (or only pieces thereof; cf. below) at 
t = 0, as discussed in jTB] for the polarized case and in [T5] for the general case. Our novel 
method here allows us to construct such solutions with arbitrary accuracy and it can hence be 
expected that this allows us to study the saddle properties of such solutions. Our main aim so 
far is to compute such a solution and hence to demonstrate the feasibility of our approach. A 
follow-up work will be devoted to the numerical construction and detailed analysis of relevant 
classes of such solutions. 

Motivated by the results in [IB], we choose the asymptotic data as follows 



u s J 1) x € [n,2ir], 

k (x) = < 1 ._i/ a! ._i/^_^ ._^, n , P**(x) = l/2, 



(x) = 0, Q**(x) 



0, x e [it, 2tt], 

g-l/Xg-l/O-a^ xe (0,7r) ; 



K{x) = k z {x) 1 A„(o:)=2. 

With these asymptotic data, the corresponding solution has a smooth Cauchy horizon at (£, x) € 
{0} x (it. 2tt) (namely where k = 1), and a curvature singularity at (i, x) € {0} x (0,7r) (namely 
where < k < 1). Note that the function k is smooth everywhere (but not analytic). Our 
analysis in Section |4.3| shows that we are allowed to set k = 1 at some points since d x Q* = 0. 
This motivates our choice of Q* ■ With this, our choice of Q** implies that the solution is polarized 
on the "domain of dependence"]^] of the "initial data" interval (w, 2tt). All data were chosen as 
simple as possible to be consistent with the constraints. 

First we repeated the same error analysis as for the previous Gowdy case, see Figure [4] For 
all the runs in the plots, we choose N = 500, At = 0.005 which guarantees that discretization 
errors are negligible in the plot. We find that our numerical method allows us to compute the 
Gowdy solution very accurately. Here, we solve the full system for (P,Q,A). 

In Figure [5] we show the numerical solution obtained from N — 1000, At = 0.0025 and 
To = —18. We plot the Kretschmann scalar at two times r = —10 and r = 0. Hence, near the 
time t = (corresponding to r = — oo), the Kretschmann scalar is large on the spatial interval 
(0, 7r) while it stays bounded at (7T, 2tt). At the later time, the curvature becomes smaller as 
expected. We also plot the remainders and of P and Q, respectively. It is instructive 
to study how the polarized region inside (tt, 2ir) gets "displaced" by the non-polarized solution. 



6 Concluding remarks 

This paper presented a new approach to the singular initial value problem for second-order 
Fuchsian-type equations. Our original motivation was to find a reliable and accurate numerical 

5 The notion of "domain of dependence" for the singular initial value problem follows from the energy estimate. 
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Figure 4: Convergence of numerical solutions corresponding to asymptotic data in Section [5. 3| 

scheme which, following [1., allowed us to evolve data numerically from the singularity. It turned 
out that the classical Fuchsian theory -despite its major successes otherwise- was not directly 
applicable to our puspose. This led us, on one hand, to develop a new approximation scheme, 
which is particularly natural to handle hyperbolic Fuchsian equations, and, on the other hand, 
to revisit the classical theory and deduce a direct existence proof for the singular initial value 
problem. Our scheme yields a reliable and accurate numerical method, referred to here as 
the Fuchsian numerical algorithm. Importantly, we demonstrated that our method applies to 
Gowdy-symmetric solutions to Einstein's field equations. 

Our method should allow us to contribute to the understanding of strong gravitational fields. 
In this direction, a particularly interesting and outstanding problem in general relativity is the 
strong cosmic censorship conjecture. Our approach allows to numerically construct, in particular, 
exceptional spacetimes, for example solutions with Cauchy horizons. This is of interest for two 
reasons, at least. First, we can learn more about the "solution space" of Einstein's field equations 
and, consequently, about the validity of general relativity as a physical theory, due to the unusual 
and sometimes physically "undesired" properties of these exceptional solutions. We are currently 



investigating the geometries of the solutions obtained in Section 5.3 in greater detail [3]. 

Second, our method allows one to study the stability of such exceptional solutions, which 
was checked here, as a first step, by perturbing the induced data on a spacelike hypersurface 
and computing the corresponding solution via the standard (backward) approach. This allows a 
systematic numerical study of the strong cosmic censorship conjecture. Note, however, that the 
strong cosmic censorship conjecture is known for Gowdy-symmetric solutions (3^1 HO]- Hence, 
in order to obtain new interesting results about this conjecture we need to apply our theory to 
more general classes of solutions; see [PJ. 

Our Fuchsian heuristics enables us to distinguish between "dominant" and "negligible" terms 
in the equations as one approaches the "singularity" . Sometimes these terms can be interpreted 
as physically interesting quantities like "kinetic" or "potential energy" . Importantly, we discov- 
ered in the present paper that the principal part of the partial different system need not, by itself, 
determines the singular behavior of the solutions. Instead, nonlinear terms (classically treated as 
lower-order source-terms) often also play an important role. This is so for Gowdy-symmetric so- 
lutions, but it is even more important for general solutions where mixmaster behavior is expected 
according to the BKL conjecture. According to this conjecture and more recent investigations 
|34l 135] , spatial derivative terms are expected to be insignificant except for exceptional points 
where spikes occur. 
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Figure 5: Numerical solutions corresponding to asymptotic data in Section pT3"| 
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